Prediction of global ionospheric VTEC maps using an adaptive autoregressive model

In this contribution, an adaptive autoregressive model is proposed and developed to predict global ionospheric vertical total electron content maps (VTEC). Specifically, the spherical harmonic (SH) coefficients are predicted based on the autoregressive model, and the order of the autoregressive model is determined adaptively using the F-test method. To test our method, final CODE and IGS global ionospheric map (GIM) products, as well as altimeter TEC data during low and mid-to-high solar activity period collected by JASON, are used to evaluate the precision of our forecasting products. Results indicate that the predicted products derived from the model proposed in this paper have good consistency with the final GIMs in low solar activity, where the annual mean of the root-mean-square value is approximately 1.5 TECU. However, the performance of predicted vertical TEC in periods of mid-to-high solar activity has less accuracy than that during low solar activity periods, especially in the equatorial ionization anomaly region and the Southern Hemisphere. Additionally, in comparison with forecasting products, the final IGS GIMs have the best consistency with altimeter TEC data. Future work is needed to investigate the performance of forecasting products using the proposed method in an operational environment, rather than using the SH coefficients from the final CODE products, to understand the real-time applicability of the method. © The Author(s) 2018. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Open Access *Correspondence: [email protected] 1 Collaborative Innovation Center for Geospatial Information Technology, No. 129 Luoyu Road, Wuhan 430079, China Full list of author information is available at the end of the article Introduction The ionosphere plays an important role in the dynamics of space weather of solar–terrestrial space. The ionosphere is important in matters of national defense, aviation security, economic development, and human life. Ionosphere monitoring using dual-frequency measurements from the Global Navigation Satellite System (GNSS) has been a topic for several decades (Komjathy 1997; De Franceschi and Zolesi 1998; Mannucci et al. 1998; Hernández-Pajares et al. 1999; Schaer 1999). GNSS provides an opportunity for long-term monitoring of the ionosphere with high accuracy and temporal and spatial resolution at relatively low cost, either in a regional context or on a global scale. The Ionosphere Associate Analysis Centers (IAACs) of the International GNSS Service (IGS) (Dow et al. 2009) have been providing reliable global ionospheric maps (GIMs) since 1998 (Hernández-Pajares et al. 2009). The IGS final vertical total electron content (TEC) maps are used for the scientific analysis of the ionosphere and practical applications. However, the IGS final GIM product is released with a time delay of approximately 2 weeks, limiting their application in real-time scenarios, including real-time precise positioning (Shi et al. 2012; Rovira-Garcia et al. 2015) and space missions, such as the Soil Moisture and Ocean Salinity (SMOS) from the European Space Agency (ESA) (Silvestrin et al. 2001; García-Rigo et al. 2011). Rapid GIM products, e.g., UQRG (ftp://newg1.upc.es/upc_ ionex) and WHUD (ftp://pub.ionosphere.cn) provided by the Technical University of Catalonia (UPC) and Wuhan University (WHU), respectively, are available with oneday latency. Meanwhile, real-time GIMs produced by IAACs will be available in the near future. However, the accuracy of real-time ionospheric products on a global scale might be limited by data availability, as the public real-time data stream is currently more concentrated in certain regions, i.e., North America, Europe, and Australia. In addition, the applications might be limited by a latency of a few seconds needed to get the real-time Page 2 of 14 Wang et al. Earth, Planets and Space (2018) 70:18 products of ionosphere. Thus, short-term predictions of global ionospheric vertical TEC (VTEC) maps are important for technological applications. Since real-time satellite orbits and clocks are available, the limiting factor in high-accuracy positioning is the ionospheric delay (Rovira-Garcia et al. 2015). Short-term predictions could be used to generate real-time global VTEC maps (Orús Pérez et al. 2010). Along with many other applications, such as automobiles, road mapping, and location-based services, short-term predictions could be used to achieve sub-meter accuracy for mass-market single-frequency receivers (García-Rigo et al. 2011). To meet the needs presented by the study of ionospheric physics and application in GNSS positioning, a few ionospheric models have been constructed, i.e., the Klobuchar model (Klobuchar 1987), the International Reference Ionosphere (IRI) model (Rawer et al. 1978; Bilitza and Reinisch 2008; Bilitza et al. 2011), and the NeQuick model (Radicella and Leitinger 2001; Nava et al. 2008). Many scholars have investigated the accuracy of these models for different regions of the world during periods of different solar and geomagnetic activities (Abdu et al. 1996; Araujo-Pradere et al. 2003; Bertoni et al. 2006; Lee and Reinisch 2006; Bhuyan and Borah 2007; Mosert et al. 2007; Adewale et al. 2011; Nigussie et al. 2013; Okoh et al. 2013; Olwendo et al. 2013; Wichaipanich et al. 2013; Wang 2016). The annual mean of the root-mean-square (RMS) of the differences between IGS GIMs and IRI predictions in 2014 was approximately 10 total electron content units (TECU, 1016 el/m2) (Wang et al. 2016). Thus, these empirical models are suitable for use in the scientific study of ionosphere behavior, which can provide predictions of the ionosphere, but they might not be appropriate for other applications that require high accuracy. Additionally, other models are built to represent the majority of the variations and the temporal–spatial distribution of the global or regional ionospheric TEC. For instance, global models are constructed by using empirical orthogonal function analysis to reproduce the major variations in TEC and the ionospheric climatology (Ercha et al. 2012; Wan et al. 2012; Mukhtarov et al. 2013). Also, regional models are studied over many countries and regions by using different methods to capture more details of the ionosphere (Bouya et al. 2010; Habarulema 2010; Chen et al. 2015; Fuller-Rowell et al. 2016; Huang et al. 2017). It is possible to obtain better VTEC maps by forecasting in the short term than those derived from empirical models. Several methods have been developed for ionospheric forecasting in recent years, such as the autocorrelation method (Muhtarov and Kutiev 1999), the autoregressive moving average (ARMA) method (Krankowski et al. 2005), a method based on neural networks (Tulunay et al. 2006), and an autoregressive model for predicting VTEC values (Karthik et al. 2012). However, many predictions are investigated at a certain location or over a regional area. In terms of global ionospheric VTEC maps forecasting, the Center for Orbit Determination in Europe (CODE), an IAAC, has provided predicted ionospheric products (oneand two-day-ahead VTEC maps, named C1PG and C2PG, respectively) for public access since 2008, via the FTP server of the Crustal Dynamics Data Information System (CDDIS, ftp://cddis.gsfc.nasa.gov/). Shortly afterward, the European Space Agency (ESA) and the Technical University of Catalonia (UPC) released their two-day-ahead VTEC maps through the FTP server of CDDIS, as well. Least-squares collocation (LSC) is used by CODE to extrapolate the coefficients of spherical harmonics (SH) for predicting VTEC maps (Schaer 1999). UPC VTEC forecasting is based on the discrete cosine transform (DCT) technique (García-Rigo et al. 2011). Moreover, a linear regression module is used to forecast the DCT coefficients and predict VTEC maps. Unlike the two methods above, the development of adaptive autoregressive modeling (AARM) for the prediction of global ionospheric VTEC maps will be presented in this manuscript. The first section of the manuscript is devoted to a detailed explanation of the AARM for ionospheric forecasting. The performance of AARM forecasting is investigated through a comparison between VTEC map predictions and IGS final products. Additionally, a comparison between VTEC predictions and external independent JASON data is conducted. Finally, conclusions are summarized in the last section. Basic methodology of ionospheric forecasting Technical description of AARM Autoregressive (AR) models have been widely used in several subject areas, such as economics (Cheng 1982), geophysics (Weiss et al. 2012), and climate change (Gu and Jiang 2005; Lee et al. 2016). First, the basic methodology of the autoregressive model (Hamilton 1994) is presented as follows: where [xt] is the time series; [a1, a2, . . . , ap] is the vector of unknown AR coefficients, which can be estimated by least square estimation (LSE); p is the order of the AR model; and [εt] is the zero-mean white noise. The usual strategy for one-step forecasting is performed using the estimated AR coefficients, as depicted in the following equation: (1) xt = a1xt−1 + a2xt−2 + · · · + apxt−p + εt (2) xn+1 = a1xn + a2xn−2 + · · · + apxn−p + εn+1 Page 3 of 14 Wang et al. Earth, Planets and Space (2018) 70:18 where xn+1 is the forecasting parameter; n is the total number of the observed time series; and εn+1 is the corresponding noise. The selected order determines the goodness of fit of the model, as well as the accuracy of the forecast. Some researchers have adopted multiple a priori computations of the model order, up to a predefined maximum order M, to select the order that gives the minimum fitting error (Costa and Hengstler 2011). The predefined maximum order M is usually determined experimentally. However, the order M may not be high enough to construct a model and forecast with high accuracy. On the other hand, if the order M is selected to be too high, a larger number of computations will be required. To avoid excessive computational cost, an adaptive approach for model order selection is presented. The modeling computation starts with a predefined minimum order N. Subsequently, modeling of the order N + 1 is also carried out. An F test is used to demonstrate whether there is a significant difference between the two models with different orders. The formula of F-statistics is as follows in Eq. (3): where RSS is the residual sum of squares; N is the predefined minimum order; and f is the degrees of freedom. In hypothesis testing, if the F-statistic is smaller than the critical value Fα, there is no significant difference between the two models. In this case, the lower-order N will be selected for modeling. Otherwise, one more modeling computation with increased order should be performed until the F-statistics is smaller than the critical value. The flowchart of AARM for extrapolation of SH coefficients is presented in Fig. 1. Forecasting global VTEC maps using AARM As a reference, the spherical harmonic coefficients produced by CODE, which are available at CODE’s FTP, are used to forecast global VTEC maps. To investigate the performance of forecasts during periods of different levels of solar activity, the SH coefficients from final ionospheric products are collected for forecasting VTEC maps in 2009 and 2015. CODE uses data from approximately 200 GNSS stations in the IGS network and other institutions. The VTEC is modeled in a solar–geomagnetic reference frame using a spherical harmonics expansion of up to degree and order 15. Piecewise linear functions are used for representation in the time domain. The time spacing of their vertices is 2 h, conforming with the epochs of the VTEC maps. (Schaer 1999). CODE divides all observations from a given day into 12 sessions, and each session contains 2 h of data. Therefore, there is (3) F = (RSSN − RSSN+1)/ (